﻿/*
 *  Copyright (C) 2002 Urban Science Applications, Inc. 
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 * 
 * 
 * for Petra.Gis adopted by asida@centrum.cz
 */


using System;
using System.Collections.Generic;
using System.Runtime.Serialization;
using Petra.Gis.Geometry.Support;
using Petra.Gis.ReferenceSystemServices;


namespace Petra.Gis.ReferenceSystemServices.Transformations
{

  
  /// <summary>
  ///		Implements the Albers projection.
  /// </summary>
  /// <remarks>
  /// 	<para>Implements the Albers projection. The Albers projection is most commonly
  /// 	used to project the United States of America. It gives the northern
  /// 	border with Canada a curved appearance.</para>
  /// 	
  ///		<para>The <a href="http://www.geog.mcgill.ca/courses/geo201/mapproj/naaeana.gif">Albers Equal Area</a>
  ///		projection has the property that the area bounded
  ///		by any pair of parallels and meridians is exactly reproduced between the 
  ///		image of those parallels and meridians in the projected domain, that is,
  ///		the projection preserves the correct area of the earth though distorts
  ///		direction, distance and shape somewhat.</para>
  /// </remarks>
  public class AlbersProjection : ProjectionTransform
  {
    private double _falseEasting;
    private double _falseNorthing;
    private double c;		  		 //constant c 
    private double e;				 //eccentricity
    private double rh;				 //heigth above elipsoid  
    private double ns0;				 //ratio between meridians
    private double lon_center; //center longitude   
    private double es = 0;

    private double const1;

    private double _semiMajor;

    private double rh1;			  // height above ellipsoid
    private double qs;			  // small q
    private double con;			  // temporary variable
    private double theta;		  // angle

    private long flag = 0;		// error flag


    private double sin_phi, cos_phi;		/* sine and cos values		*/


    #region constructors



    public AlbersProjection()
      : base()
    { }




    /// <summary>
    /// Creates an instance of an Albers projection object.
    /// </summary>
    /// <param name="parameters">List of parameters to initialize the projection.</param>
    /// <remarks>
    /// <para>The parameters this projection expects are listed below.</para>
    /// <list type="table">
    /// <listheader><term>Items</term><description>Descriptions</description></listheader>		<item><term>latitude_of_false_origin</term><description>The latitude of the point which is not the natural origin and at which grid coordinate values false easting and false northing are defined.</description></item>
    /// <item><term>longitude_of_false_origin</term><description>The longitude of the point which is not the natural origin and at which grid coordinate values false easting and false northing are defined.</description></item>
    /// <item><term>latitude_of_1st_standard_parallel</term><description>For a conic projection with two standard parallels, this is the latitude of intersection of the cone with the ellipsoid that is nearest the pole.  Scale is true along this parallel.</description></item>
    /// <item><term>latitude_of_2nd_standard_parallel</term><description>For a conic projection with two standard parallels, this is the latitude of intersection of the cone with the ellipsoid that is furthest from the pole.  Scale is true along this parallel.</description></item>
    /// <item><term>easting_at_false_origin</term><description>The easting value assigned to the false origin.</description></item>
    /// <item><term>northing_at_false_origin</term><description>The northing value assigned to the false origin.</description></item>
    /// </list>
    /// </remarks>
    public AlbersProjection(ProjectedCoordinateSystem proj)
      : this(proj, false) { }





    public AlbersProjection(ProjectedCoordinateSystem proj, bool isInverse)
      : base(proj, isInverse) { }

    #endregion constructors



    protected override bool parametersValid()
    {
      return Parameters.ContainsKey("latitude_of_origin")
        && Parameters.ContainsKey("central_meridian")
        && Parameters.ContainsKey("standard_parallel_1")
        && Parameters.ContainsKey("standard_parallel_2");
    }
  
  
  
    protected override void updateTransform()
    {
      double sin_po, cos_po;	// sin and cos values
      double ms1;				// small m 1
      double ms2;				// small m 2
      double qs0;				// small q 0
      double qs1;				// small q 1
      double qs2;				// small q 2

      double lat0 = Parameters["latitude_of_origin"] * AngularParameterUnit.RadiansPerUnit;
      double lon0 = Parameters["central_meridian"] * AngularParameterUnit.RadiansPerUnit;
      double lat1 = Parameters["standard_parallel_1"] * AngularParameterUnit.RadiansPerUnit;
      double lat2 = Parameters["standard_parallel_2"] * AngularParameterUnit.RadiansPerUnit;

      if (!Parameters.ContainsKey("false_easting"))
        _falseEasting = 0;
      else
        _falseEasting = Parameters["false_easting"] * DistanceParameterUnit.MetersPerUnit;

      if (!Parameters.ContainsKey("false_northing"))
        _falseNorthing = 0;
      else
        _falseNorthing = Parameters["false_northing"] * DistanceParameterUnit.MetersPerUnit;

      lon_center = lon0;

      if (Math.Abs(lat1 + lat2) < ProjMath.EPSLN)
        throw new ArgumentException("Equal latitudes for St. Parallels on opposite sides of equator.");

      e = Projection.GeographicCoordinateSystem.Datum.Ellipsoid.Eccentricity;
      es = Projection.GeographicCoordinateSystem.Datum.Ellipsoid.e2;

      _semiMajor = Projection.GeographicCoordinateSystem.Datum.Ellipsoid.SemiMajorAxis
        * Projection.GeographicCoordinateSystem.Datum.Ellipsoid.AxisUnit.MetersPerUnit;

      const1 = 1 - .5 * (1.0 - es) * Math.Log((1.0 - e) / (1.0 + e)) / e;

      ProjMath.Sincos(lat1, out sin_po, out cos_po);
      con = sin_po;

      ms1 = ProjMath.Msfnz(e, sin_po, cos_po);
      qs1 = ProjMath.Qsfnz(e, sin_po, cos_po);

      ProjMath.Sincos(lat2, out sin_po, out cos_po);

      ms2 = ProjMath.Msfnz(e, sin_po, cos_po);
      qs2 = ProjMath.Qsfnz(e, sin_po, cos_po);

      ProjMath.Sincos(lat0, out sin_po, out cos_po);

      qs0 = ProjMath.Qsfnz(e, sin_po, cos_po);

      if (Math.Abs(lat1 - lat2) > ProjMath.EPSLN)
        ns0 = (ms1 * ms1 - ms2 * ms2) / (qs2 - qs1);
      else
        ns0 = con;
      c = ms1 * ms1 + ns0 * qs1;
      rh = this._semiMajor * Math.Sqrt(c - ns0 * qs0) / ns0;

    }

    #region methods




    /// <summary>
    /// Converts coordinates from lon, lat to x,y
    /// </summary>
    /// <param name="p">coordinate to transform</param>
    protected override void geographicToProjectedN(ICoordinate p)
    {

      //// get to radians
      //p.lon *= _geogUnit.RadiansPerUnit;
      //p.lat *= _geogUnit.RadiansPerUnit;

      ProjMath.Sincos(p.lat, out sin_phi, out cos_phi);
      qs = ProjMath.Qsfnz(e, sin_phi, cos_phi);
      rh1 = _semiMajor * Math.Sqrt(c - ns0 * qs) / ns0;
      theta = ns0 * ProjMath.AdjustLongitude(p.lon - lon_center);

      p.X = rh1 * Math.Sin(theta) + _falseEasting;
      p.Y = rh - rh1 * Math.Cos(theta) + _falseNorthing;

      //// convert to output units
      //p.X /= _projUnit.MetersPerUnit;
      //p.Y /= _projUnit.MetersPerUnit;

    }





    /// <summary>
    /// transform coordinate in projected units x,y to angular lon,lat
    /// </summary>
    /// <param name="p">coord to transform</param>
    protected override void projectedToGeographicN(ICoordinate p)
    {

      ////convert source coordinate to meters
      //p.X *= _projUnit.MetersPerUnit;
      //p.Y *= _projUnit.MetersPerUnit;

      p.X -= _falseEasting;
      p.Y = rh - p.Y + _falseNorthing; ;
      if (ns0 >= 0)
      {
        rh1 = Math.Sqrt(p.X * p.X + p.Y * p.Y);
        con = 1.0;
      }
      else
      {
        rh1 = -Math.Sqrt(p.X * p.X + p.Y * p.Y);
        con = -1.0;
      }
      theta = 0.0;
      if (rh1 != 0.0)
        theta = Math.Atan2(con * p.X, con * p.X);
      con = rh1 * ns0 / this._semiMajor;
      qs = (c - con * con) / ns0;
      if (e >= ProjMath.EPSLN)
      {
        //con = 1 - .5 * (1.0 - es) * Math.Log((1.0 - e) / (1.0 + e)) / e;
        if (Math.Abs(Math.Abs(const1) - Math.Abs(qs)) > .0000000001)
        {
          p.lat = ProjMath.Phi1z(e, qs, out flag);
          if (flag != 0)
            throw new ApplicationException();

        }
        else
        {
          if (qs >= 0)
            p.lon = ProjMath.HALF_PI;
          else
            p.lat = -ProjMath.HALF_PI;
        }
      }
      else
      {
        p.lat = ProjMath.Phi1z(e, qs, out flag);
        if (flag != 0)
          throw new ApplicationException();

      }

      p.lon = ProjMath.AdjustLongitude(theta / ns0 + lon_center);

      //// update to target units
      //p.lon /= _geogUnit.RadiansPerUnit;
      //p.lat /= _geogUnit.RadiansPerUnit;

    }





    /// <summary>
    /// Returns the inverse of this projection.
    /// </summary>
    /// <returns>reverse of the current projection.</returns>
    public override SRSTransform Inverse()
    {
      return new AlbersProjection(Projection, !IsInverse);
    }





    public override Dimension TransformDimension
    {
      get { return Dimension.dim2D; }
    }


    #endregion methods




    public static string[] ProjectionNames
    {
      get { return new string[] { "alberts", "alberts projection", "alberts_projection" }; }
    }


  
  }

}
